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This paper presents a semi-analytic method for computing frequency de- 
pendent means, variances, and failure probabilities for arbitrarily large- 
order closed-loop dynamical systems possessing a single uncertain parame- 
ter or with multiple highly correlated uncertain parameters. The approach 
will be shown to not suffer from the same computational challenges associ- 
ated with computing failure probabilities using conventional FORM/SORM 
techniques. The approach is demonstrated by computing the probabilistic 
frequency domain performance of an optimal feed- forward disturbance re- 
jection scheme. 


N omenclat ure 

C = cubic spline approximation of F (when U -space is one dimensional) 

Ci = cubic polynomial that represents C on interpolation interval number i 

F = 17- space counterpart of /: F(U) = f(T~ l (U)) so Y = F(U) 

Fx = cumulative distribution function of X 

Fy — cumulative distribution function of V 

f = function relating X and Y:Y — f(X) 

G — counterpart of g on {/-space: G(u) = g(T~ 1 (u)) 

g — limit state function on A-space 

Pf = probability of failure 

T = transformation from A r -space into {/-space 

U — scalar or vector of uncorrelated standard normal random variables 
u* = most probable point (MPP) of failure 

X = scalar or vector of random variables 

x = generic point in A'-space 

Y = scalar or vector of random variables dependent on A : Y = /(A) 

Y a = approximation of Y: Y a ~ C(U) = C(T( A)) 
y = generic point in Y -space 

0 — ||u*||, the reliability index 

$ =1 dimensional standard normal cumulative distribution function 
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I. Introduction 

Stability and performance robustness are two of the fundamental metrics 
that must be considered in control system analysis and design. The first sys- 
tematic procedures for analyzing the robustness of closed-loop dynamical sys- 
tems were developed by H. W. Bode 1 and H. Nyquist 2 . These researchers were 
instrumental in developing some of the fundamental concepts for the analysis 
and design of feedback control systems and their frequency domain methodolo- 
gies for single-input/single-output (SISO) systems that are still in use. Today 
in the digital world with high performance computers and powerful numeri- 
cal software, the SISO concepts for analyzing robustness have been generalized 
to include multi- input /multi-output (MIMO) systems using loop-gain singular 
value analysis 3 , structured singular- values 4,5 , and //-analysis 6 . Unfortunately, 
these conventional approaches in many cases result in very conservative esti- 
mates of stability and performance and rely upon a mathematical definition of 
the uncertainties that may not necessarily reflect the true nature of the param- 
eter and model-form uncertainties. In an attempt to overcome these limitations, 
researchers are now investigating methods for analyzing the stability of feedback 
systems with parameter uncertainties defined in terms of probability distribution 
functions. These approaches have been referred to as “probability of stability” , 
see Refs. 7 9. Initial work in this area relied upon Monte Carlo-based methods 
to assess the effect of probabilistic parametric uncertainties on system stability. 
Although viable in many cases, simulation based approaches have computational 
weaknesses when analyzing very low probability events, such as those related 
to aerospace vehicle stability.* 1 More recently, researchers 10,11 have explored the 
use of first-order and second-order reliability methods (FORM and SORM) for 
computing probabilistic performance robustness in terms of root mean square 
(RMS) output response and probabilistic stability analysis in terms of the sys- 
tem’s eigenvalues. Although very important in the analysis of performance ro- 
bustness of controlled systems, RMS response is an integrated quantity and does 
not provide the control system designer with important response information 
about the spectral nature of dynamical systems. 

This paper extends previous probabilistic works by considering the full spec- 
tral domain of a closed-loop dynamical system and demonstrates the com- 
putation of means, variances, and probabilistic confidence intervals across a 
broad range of frequencies. Furthermore, this paper will highlight some of the 


a See also the paper “Multicriteria Optimization Methods for Design of Aircraft Control Systems” 
by Albert A. Schy and Daniel P. Giesy in the book Multicriteria Optimization in Engineering 
and in the Sciences (Plenum Press, New York, 1988; pp. 225-262) edited by Wolfram Stadler, 
and references [16] and [17] from that paper. In the example Schy and Giesy present, the control 
effectiveness and sideslip aerodynamic derivatives of the Space Shuttle are the uncertain parameters. 
These appear in the equations of lateral motion which are, in turn, used to model a lateral stability 
augmentation system (SAS) based on state feedback. Using the feedback gains as design variables, 
the SAS is designed by simultaneously optimizing multiple criteria. The multiple criteria are based 
on the eigenstructure of the closed loop lateral equations of motion and on the roll response, peak 
sideslip, peak control deflection, and peak control rate values over the 6 second roll maneuver 
duration. In order to make the design process Stochastic Insensitive , not only are the nominal 
values of the multiple criteria optimized but also estimates (based on the assumption that the design 
criteria are linear functions of the uncertain parameters) of the probabities that each criterion will 
fall within designer specified desirable bounds. 
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computational challenges associated with computing failure probabilities using 
conventional FORM/SORM techniques. Additionally, this paper will present a 
semi- analytic method for computing means, variances, and failure probabilities 
for arbitrarily large-order systems possessing a single uncertain parameter or 
with multiple highly correlated uncertain parameters. The methodology will be 
demonstrated by computing the probabilistic performance of an optimal feed- 
forward disturbance rejection scheme with a single uncertain parameter. The 
case of multiple highly correlated uncertain parameters is dealt with by approx- 
imating the uncertain parameters as functions of a single uncertain parameter 
and applying the techniques developed in this paper. 

II. Mathematical model of uncertainty 

In a standard mathematical formulation used for uncertainty modeling 12 , a 
scalar or vector A" of random variables is used to represent one or more uncertain 
parameters. A real valued function g , called the limit state function, divides the 
uncertain domain into success (#(X) > 0) and failure {g{X) < 0) regions. The 
boundary between these regions, defined by g(X) = 0, is called the limit state 
surface. 

One statistic of interest in design and analysis problems is the probability 
of failure, Pf — P[#(X) < 0]. This is given exactly by 

P/=/ 1 dF x (1) 

Jg(X)<0 

where Fx is the cumulative distribution function (CDF) of A". However, eval- 
uation of this integral is usually computationally intractable and so, histori- 
cally, recourse has been made to approximation techniques such as Monte Carlo, 
FORM, or SORM. 

To approximate failure probability using FORM or SORM, X is transformed 
by a transformation T so that U = T(X) is a vector of uncorrelated standard 
normal random variables. The X-space limit state function g is transformed 
to a [/-space limit state function G by G(m) = g(T~ l (u)). The closest point 
to the origin on the [/-space limit state surface G(u) = 0, called the most 
probable point (MPP) of failure and denoted u*, is determined by solving the 
optimization problem “minimize ||tx[| subject to the constraint G(u) = 0”. The 
norm of u* is denoted by /? and called the reliability index. FORM approximates 
Pf by replacing the [/-space limit state surface G(u) = 0 by the hyperplane 
tangent to that surface at u*. Then Pf is approximated by 1 — $(/?) where 
is the CDF of the one dimensional standard normal random variable. This Pf 
approximation is exactly the probability that U is in the half-space determined 
by the tangent hyper plane that does not contain the origin. This is appropriate 
if G(0) > 0. If G(0) < 0, Pj is approximated by <3>(/3), while if G(0) = 0, the 
FORM approximation to Pf is 0.5. 

The SORM approximation makes a second order approximation to the limit 
state surface at the MPP and so can be expected to provide a better approx- 
imation to the probability of failure. However, a cost is incurred both in the 
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need to have curvature information for the limit state function and in a more 
complicated evaluation of the probability of a standard normal random variable 
falling on one side of a parabolic or quadratic surface 12-14 . The present study 
made no use of SORM. 

FORM has been successfully applied in problems of real engineering interest, 
such as structural component reliability problems. These successes occur when 
the probability of failure is small and the failure set (at least in U- space) does 
not have too pathological a geometry. However, if there are multiple points on 
the limit state surface that are not in close proximity to each other and whose 
norm is nearly /3, or if u* is fairly close to 0 and the limit state surface has 
substantial curvature at u *, FORM can provide a poor approximation to the 
failure probability. 

A simple example of this problem is given in Fig. 1. In this one-dimensional 
example, the failure set consists of two semi-infinite intervals. The limit state 
“surface” is zero-dimensional, i.e., a collection of (in this case 2) discrete points; 
namely -1.0 and +1.1. Thus the MPP is the closer of these two to 0, namely -1.0. 
The FORM estimation of the failure probability is based on the probability P x 
that U falls in the left component of the failure set and ignores the contribution 
P 2 of the right component of the failure set. It thus misses the mark by about 
46%. 


III. Cumulative Distribution Functions of Dependent 

Random Parameters 

Mathematically, an engineering system may be viewed as a vector-valued 
function of several variables. The independent variables could include design 
parameters such as component dimensions and controller gains and physical 
givens such as material properties. The dependent variables could include a wide 
variety of metrics of interest such as total mass, load capacity, stability mar- 
gins, etc. If any of the independent variables are uncertain, arid are therefore 
represented by random variables, the dependent variables also become uncer- 
tain. It is of interest in analyzing the system to determine what is the CDF of 
each dependent variable. Once this CDF is known, it can be used to compute 
statistics of interest such as mean, variance, confidence intervals, etc. This CDF 
is determined by the function that describes the system and the joint CDF of 
the independent variables. The next topic of interest is determining the CDF of 
a dependent variable; multiple dependent variables would be treated one at a 
time. 

Consider a random parameter Y that is a function of other random param- 
eters: Y = f{X) . The probabilistic behavior of Y is described by its CDF, Fy , 
which is determined by f{X) and the CDF Fx of X: 

Fy(v) = P{Y <y\ = P[f(X) < y\ (2) 

This can be modeled in the same mathematical pattern as is used in the FORM 
approximation. If the limit state function g is defined as g(x) = f(x) — y, the 
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Fig. 1. A Simple Limit State Function with Multiple MPPs 

“probability of failure” computed using this g approximates Fy(y). (It is proba- 
bility a misnomer to continue to use the historical phrase “probability of failure” 
to describe the results of this calculation. The event [Y < y] is a “failure” only 
in that the random variable Y fails to exceed the sample value y. Notwithstand- 
ing, we will not introduce new terminology to use when traditional uncertainty 
analysis tools are used to approximate CDFs instead of failure probabilities.) 
This probability may now be approximated by FORM or SORM. Note that for 
each y for which the value of Fy(y) is to be calculated, an optimization prob- 
lem must be solved to find the MPP for the corresponding limit state function. 
The computational intensity of this task is mitigated by using the MPP from 
one value of y as a starting point in the calculation of the next - although a 
small change in y might result in a large jump in the location of the MPP, it 
is reasonable to expect that in most cases the MPP will change smoothly with 
changing y, and using the MPP from a neighboring problem will provide a good 
approximation to the MPP of the present problem leading to fast convergence 
of the optimization process. 
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One CDF of interest in the present study was the CDF of the Bode mag- 
nitude plot of a system with at least one random parameter. This adds a new 
dimension of complication - the CDF is really a family of functions, one for each 
frequency. If the FORM technique is to be used to approximate these CDFs. it 
would take advantage of information from calculations made at one frequency 
value to give a head start to optimization problems being solved at a near by 
frequency. 

However, when this technique was applied in practice, validity checking with 
Monte Carlo simulation revealed that some of the problems that can beset 
FORM approximation were actually occurring. Failure probability sets of more 
than one component, such as is depicted in Fig. 1 , were occurring in the Bode 
magnitude data considered in section V. This resulted in FORM based approx- 
imations to the CDF which were unacceptably inaccurate. When the uncertain- 
ties could be expressed in terms of a single random variable, an alternative to 
FORM was found to approximate the CDFs. 

IV. The one-dimensional case 

The case considered here is the one in which X and therefore U are scalar 
random parameters. Recall that U has the standard normal distribution. We 
wish to find the CDF of Y = f(X). In the application of interest in this paper, 
f(X) will be the value of the Bode magnitude plot at a fixed value of frequency 
of some system in which X is a parameter. b When X is transformed into U- 
space, Y becomes a function of U: Y = F(U) where F(u) = f(T~ 1 (u)). For each 
y, Fy = P[Y <y] = P[F(U) < y\. Suppose that the function F is reasonably 
well-behaved (as could be expected in engineering applications). Then, if the set 
{u\F(u) < y} is not empty, it is a collection of intervals (possibly some finite, up 
to two semi-infinite, or just the whole real line). The boundary points of these 
intervals come from the solution set of the equation F(u) = y. If a complete 
set of solutions to this equation can be found, and if F is nice enough that the 
solution set contains only finitely many points, then these points divide the real 
line into a finite number of intervals on the interior of each of which we have 
that F is either strictly greater than y or strictly less than y. If U\ < U2 are 
consecutive solutions of F(u) = y , then it can be determined whether F(u) > y 
or F(u) < y on the open interval (1x1,112) by evaluating F at a single interior 
point or evaluating F' at either endpoint. Adding up the probabilities that U 
lies in each interval of the set of intervals where F(u) < y produces the exact 
value of F y (y ); i.e., the integral in Eq. ( 1 ) has been evaluated. 

There are two problems with this. 

1 . For a general function F, it is difficult to be sure that a numerical pro- 
cedure has found all solutions to F(u) = y. To see the problem if the 


b The following discussion has been set in U- space. It could just as easily have been set in X- 
space. One would just use the original f instead of the transformed F. In place of using 4> to find 
the probability that U falls in a given interval, one would use F\ to find the probability that X 
falls in the corresponding X-space interval. 
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solution set is incomplete, presume that a numerical procedure has cal- 
culated a collection of solutions to F(u) = y . Presume also that u\ < ?/ 2 
are adjacent members of this collection of solutions. Presume further that 
wo is a solution of F(u) ~ y which the numerical procedure has failed 
to detect and which falls between u\ and i/ 2 . The most common behav- 
ior of F in the neighborhood of uq is that it changes signs as u passes 
through no. This means that the interval (u^u^ intersects both of the 
events [F(u) > y] and [F(u) < y] with positive probability. Thus, either 
including or excluding P[u\ < U < U 2 } in the calculation of P[F(£/) < y] 
introduces an error in that calculation. 

2. Any attempt to calculate all of the solutions of F(u) = y would require 
many (potentially expensive) evaluations of the function F. 

The solution to these problems used in this study was to use a Response Surface 
method: F is replaced by a cubic spline approximation C and Y = F(U) is 
approximated by Y a = C(U). An effective and efficient procedure is devised for 
calculating all of the solutions of C(u) — y over the support interval of C. As 
above, this provides the necessary information to calculate P[Y a < y] which is 
then taken as an approximation for P[Y < y]. (The contribution to P[Y a < y ] 
from points u outside the support of the cubic spline C is actually ignored, but 
this error can be made negligible by making the support of C big enough. For 
example, if the support of C is [-5, 5], so that the ±5<r range of U is covered, 
then the ignored probability is less than 5.8 x 10 -7 .) 

The repeated evaluations of F(u) that would be required for each value of 
y for which the equation F(u) = y was to be solved are now replaced by the 
evaluation of F at a collection of points (the so called knots of the cubic spline) 
that are chosen to produce a good approximation of F by C. Instead of doing 
many evaluations of F for each y, the calculation of F at the knots of C need be 
done only once, the coefficients of each of the cubic polynomials which represents 
C on the interpolating intervals are also calculated only once, and this defining 
information for C is reused for each value of y. 

At any value of y for which the CDF of Y is to be approximated, the con- 
tribution to the CDF of Y a = C(U) on the support of C can now be evaluated 
exactly. The solutions (if any) of C(u) = y are found on each interpolating in- 
terval. If Ci is the cubic polynomial that C uses to interpolate F over the i ^ 
sub-interval, then the cubic equation Ci(u) — y — 0 is solved for its three roots, 
and any that are real numbers actually falling within the subinterval are 
retained as solutions of C(u) = y. After some alternative ways of finding these 
roots were studied (use of the Matlab® polynomial root finder (roots), use of 
Sturm sequences to avoid actually calculating the roots of a cubic if they did not 
lie in the target interval), the most efficient way found to calculate the solutions 
of C(u) — y = 0 was to use the cubic formula 15 ” 17 to calculate the solutions of 
Ci{u) — y = 0 for every i. 

Since it is possible that individual polynomials of the cubic spline degenerate 
to lower order than 3, the polynomial solver implemented to locate the zeros of 
C(u) — y must be able to solve linear, quadratic, and cubic equations. This was 


7 

American Institute of Aeronautics and Astronautics 



done in the present study with software which provided numerical approxima- 
tions to the roots of polynomial equations of degree 3 or less which had relative 
errors on the order of a small multiple of machine epsilon. This was adequate 
solution accuracy in most cases. However, numerical problems can occur at the 
knots and at derivative zeros (including local maxima and minima) of the cubic 
spline, C, and special care must be exercised at these points. 

The problem which can occur if a solution to Ci(u) = y falls at one of the 
knots which constitute the boundary of the subinterval is that, when the 
solver has approximated the solution, a small error in its computed value might 
move it out of that subinterval. This results in the solution logic rejecting the 
value as a solution of C(u) = y. If that happens on both intervals which share 
that knot, the solution point will be overlooked. If it happens in neither interval, 
then the knot is entered twice in a list of computed solutions. The problem of 
duplicate solutions is fixed by post processing the list of solutions to remove 
duplicates. This also deals with any problem the probability calculation might 
have with any multiple zeros of C(u) = y detected by the polynomial solver. 
The problem of a solution at a knot being overlooked is solved by just directly 
evaluating the cubic spline at each knot and adding any knots Uj for which C{u t ) 
is within some small error tolerance of y to the solution list. Any knot detected 
as a solution point is deflated out of its cubic polynomial before any remaining 
solutions in that subinterval are sought. 

The problem which occurs at a point v! where C{u f ) = y and C f (u f ) = 0 
is that u' ought to appear as a multiple solution to C{u) = y. However, as a 
result of round-off error while evaluating the formulae which give these answers, 
what might come out of the solution process is 2 or 3 distinct real or complex 
approximations to the multiple solution of C(u) = y at u' which are closely 
clustered about u f . If all these approximations are real numbers, this is not a 
problem with regard to computing P[Y a < y\. The spurious intervals introduced 
by partitioning the real line using these approximate solutions are so short that 
any contribution to P[Y a < y] which comes from them is negligible. However, 
since only real roots are sought, any of these approximations to u f which are 
complex by virtue of having small, but nonzero, imaginary parts are discarded. 
To avoid this, the solution software is structured so that the user can pass in a 
collection {uj} of zeros of C' together with the values { yj = C(uj)} which C 
takes at these derivative zeros. If the software trying to solve C(u) = y detects 
that any of the yj are equal to y, the corresponding Uj is included in the solution 
set. The same software used to solve equations of the form C(u) = y can be 
used to solve C'(u ) = 0 and so locate the derivative zeros of C. 

(The authors are aware that, normally, any test involving comparing two 
floating point numbers for equality, as was just suggested in comparing the 
value y at which the CDF of Y a is being computed to members of the set {y 3 } 
of zeros of C", is questionable from the point of view of numerical analysis. 
That notwithstanding, it is done here because of experience based on a peculiar 
property of the CDF of a random variable dependent on another through a 
function which has derivative zeros. If any random variable is dependent on 
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another through a function which possesses points with horizontal tangent lines, 
then at the corresponding point in the function image space the CDF of the 
dependent random variable has a vertical tangent. The practice followed in this 
study was to refine the grid of points at which the dependent variable CDF was 
being calculated in the vicinity of such vertical asymptotes to better resolve the 
effect of this sudden rise in the CDF on such quantities as the mean and variance 
of the dependent random variable which are determined from the CDF. In this 
refinement, the point of actual vertical tangency was specifically included. At 
such points, the equality test between y and at least one of the yj is met.) 

Once this is done, all solutions of C(u) = y can be reliably computed, and 
(except for the neglected portions outside the support of C) the CDF of Y a is 
computed exactly at y . 

V. Application to Optimal Feed-forward Control 

A model of the flexible beam test article, shown in Fig. 2, will be used 
to validate the probabilistic analysis methodology developed in this work. The 
system consists of a very flexible thin aluminum blade, approximately one-meter 
long, attached at its base to a hub motor. The hub motor is the primary control 
effector for the system. At the tip of the beam is a reaction wheel that serves as 
a disturbance generator to provide harmonic imbalance forces. The test article 
has nine sensors that may be used in any combination for either feedback or 
performance output monitoring. The finite element method was used to model 
this system and employed a simple implementation of Euler-Bernoulli planar 
beam elements. 0 For a complete description of the flexible beam test article as 
well as the optimal feed- forward disturbance rejection methodology see Ref. 18. 

A. Probabilistic Analysis 

The analysis approach presented in this paper has been used to evaluate the 
probabilistic robustness of the system’s pointing performance for two different 
control system architectures. Two architectures were chosen as a representative 
example of the design choices facing control engineers and to illustrate how prob- 
abilistic performance analysis may be used to enhance the design process. The 
first control system that is considered is purely a Linear Quadratic Gaussian 
(LQG) feedback controller using the hub actuator to mitigate the effect of reac- 
tion wheel imbalance forces on the tip acceleration. The second control system 
utilizes the same LQG feedback structure as the first, but adds a feed-forward 
path to improve the ability to reject the reaction wheel imbalance forces. Al- 
though many control system performance metrics can be evaluated using the 
approach presented in this paper, these architectures were evaluated solely in 
terms of their probabilistic Bode magnitude response. 

The uncertain parameter considered in this study was a normally distributed 
mass density property, i.e., mass density of the beam elements used in the finite 


c The finite element model was used to produce a modal model using one rigid body mode 
(rotation about the hub motor axis) and the first seven flexible modes. This gave a 16 state plant 
model. With the addition of a full state estimator, the analysis model was a 32 state closed loop 
system. 
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source 


Experimental Test Article 


Fig. 2. Flexible Beam Test Article 

element model. d This multi-disciplinary problem required integrating structural 
analysis, control system analysis and probability tools. This integration was 
accomplished by developing all software tools in Mat lab®. The results of the 
probabilistic robustness analysis are shown in Figs. 3 and 4. Figure 3 presents 
the probabilistic performance for both control architectures in terms of their 
mean responses together with 99% confidence intervals. The computation of the 
mean as well as the confidence intervals requires the CDF of the Bode magnitude 
response for each frequency of interest. The results shown in this paper were 
generated using 200 frequency points with 40 points per CDF analysis and 50 
evaluations over a ±5cr range of the uncertain mass density variable used to 
generate the cubic spline interpolating functions. 

Once given the frequency dependent CDF, higher order statistics may also be 
easily computed. Fig. 4 is an example of a second-order moment (variance) of the 
system’s acceleration response. This figure indicates that the system with feed- 
forward control has (in general) smaller variances than the system with feed- 
back only. From this, we can conclude that the robustness of the feed-forward 
approach is not adversely affected by uncertainty in the mass density. More 
specific conclusions regarding the smaller variances are outside the scope of this 
paper, but credible variance information should be seen as a significant enhance- 
ment over the information resulting from convention robust control analysis and 
is a particularly valuable metric to be considered in design. 


'’Beam element mass density mean was 2850 Kg/m 3 . Variance is defined by the requirement that 
3 a = 285 Kg/m 3 , which is 10% of the mean. 
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Probabilistic Closed-Loop Performance 



Fig. 3. Probabilistic Analysis Results 


The results in Fig. 3 show that for the frequency range between 5-45 Hz the 
control system with feed-forward has substantially better mean value perfor- 
mance than the system with feedback only. This assessment is true not only in 
terms of the mean value performance, but also for all possible responses within 
the 99% confidence intervals. In the frequency ranges above and below this 
range, the probabilistic analysis shows that the two architectures have overlap- 
ping performances. This information is of great value for considering resource 
allocation trade studies. For example, if the system’s dominant operating range 
is in the 50-70 Hz range, then the more expensive feed-forward architecture pro- 
vides little or no performance advantage over the standard LQG controller. The 
real value of the approach is that the quantitative behavior of different control 
system architectures can be accessed in the presence of uncertainties that are 

• truly consistent with the physical phenomenon. 

As a final note, the overall central processor unit time required to compute 
the data presented in Figs. 3 and 4 was 113 seconds for each control architecture. 

• Timing data was obtained on an AMD 2600XP processor system. 
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Tip Acceleration Variance, m/s? 



Frequency, Hz 


Fig. 4. Acceleration Output Variance 
B. Monte Carlo Verification 

Since some approximation has been involved in the preceding analysis (replac- 
ing, at each frequency, the actual function by which the Bode magnitude re- 
sponse is determined by the uncertain mass density with a cubic spline approxi- 
mation; ignoring the chance that the mass density falls outside the ±5cr range), 
a Monte Carlo study was performed to validate the results for the controller 
which included feed- forward. A set of 60,000 samples of the uncertain mass 
density was generated and the Bode magnitude analysis was performed at each 
of these 60,000 sample values at each of the 200 frequency points considered 
previously. For each frequency point, the mean and variance of the Monte Carlo 
results were calculated. 

The Monte Carlo results agree so well with the analytical results using the 
cubic spline approximation that, when they were plotted together, the differ- 
ence between the two is all but lost in the limits of the resolution of the plotting * 

device. The actual error curves are shown in Figs. 5 and 6. The discrepancy 
between the mean as calculated by the cubic spline approximation and as cal- 
culated by Monte Carlo is never more than 0.2% of the maximum magnitude 
Monte Carlo value. The discrepancy between variances is never more than 1.6%. 
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Error in Tip Acceleration Mean Value, dB 



Fig. 5. Difference between Monte Carlo and Semi-analytic - Mean 


This validates the technique presented here. 

The Monte Carlo calculation required 40 hours of computer time. Compar- 
ing this computation time with the 113 seconds required for the cubic spline 
approximation method shows the advantage in efficiency of the technique pre- 
sented here. 


VI. Conclusions 

This paper has presented an approach for computing frequency dependent 
cumulative distribution functions for closed-loop dynamical systems. The ap- 
proach has been shown to not suffer from the same computational challenges as- 
sociated with computing failure probabilities using conventional FORM/SORM 
techniques. The computation of means, variances, and probabilistic confidence 
intervals has been demonstrated on a multi-disciplinary model of a laboratory 
test article for a broad range of frequencies. The advantages of this approach 
over Monte Carlo methods have been demonstrated. 
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